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Abstract 

In the study of time-dependent waves, it is computationally expensive 
to solve a problem in which high frequencies (shortwaves, with wavenum- 
ber k = fcmax) and low frequencies (longwaves, near k = fcmin) mix. Con- 
sider a problem in which low frequencies scatter off a sharp impurity. The 
impurity generates high frequencies which propagate and spread through- 
out the computational domain, while the domain must be large enough 
to contain several longwaves. Conventional spectral methods have com- 
putational cost proportional to 0(fcinax/fcmin log(fcmax/fcmin)). 

We present here a multiscale algorithm (implemented for the Schrodinger 
equation, but generally applicable) which solves the problem with cost 
(in space and time) 0(fcmaxi log(fcnmx/fcmin) log(fcmaxi)). Here, L is the 
width of the region in which the algorithm resolves all frequencies, and is 
independent of fcmin. 



1 Introduction and Definitions 

Consider the time dependent Schrodinger equation with {x,t) e M^+^: 

idt^Pix, t) = h(l/2)A + V{x)] i;{x, t) (1.1a) 

^{x,t = Q) = Mx) (1.1b) 

In this work we restrict ourselves to the case TV = 1, although there is nothing 
intrinsic to our method which requires this. 

Suppose that for the chosen initial data, ^p{k,t) remains localized between 
the frequencies fc,„i„ and fcmax (with k^in <C fcmax), apart from some small error. 

To solve this numerically, one typically truncates the domain to a finite 
region (the interaction region) and imposes some sort of open boundary con- 
ditions. The natural question to ask at this point is "what is the interaction 
region?" 

Assume V{x) is supported near x — 0. Waves have an interaction length 
proportional to their wavelength, implying that waves with wavelength A = 
27r/fc interact with V{x) over a distance 0{\) = 0{k^^). Thus, if frequencies 



are bounded below by /cmin, the interaction region has width at least 0(fc~?j^). 
Accurate resolution of high frequency waves requires at least two samples per 
period of the highest frequency fcmax, or a sampling rate O(fcmax)- The number 
of samples is ©(femax/^min) (sampling rate times box width). In phase space 
terms, the solution is being computed on the region 

\{x,k):\x\<-^,\k\<kra^A, (1.2) 

and the number of samples required is proportional to the phase space volume 
of this region. 

However, high frequency waves (with, e.g. k = fcmax/2) do not interact with 
the potential past x = 0{k~^). So the interaction region in phase space is: 

X = {{x, k) : < C/ \x\ , \k\ < k^^}. (1.3) 

Assuming ipik,t) contains little mass in [— fcmin, fcmin], we can truncate at x = 
0{k~lJ. This tnmcated region has volume 0{kynax^og{k^l^)) rather than 
0(fcmaxfcj„i„), which is asymptotically smaller (the constant of proportionality 
is proportional to C^/^ and has units of length). This observation suggests that 
ordinary spectral methods are inefficient for studying interactions, and that 
more efficient numerical methods are possible. 

In this work we present an algorithm exploiting this with spectral accu- 
racy. The computational complexity is 0{Lkmax log(fcmax/A;min) log(iA;inax)) per 
timestep, rather than 0(fci„ax/fcmin) for normal spectral methods. Additionally, 
being a Fourier-based spectral method, it can be easily modified to treat other 
wave equations. Here, L is chosen so that the high frequency components of T 
are entirely contained [—L,L] (i.e. if (x, fcmax/2) S I, then x G [—L,L]). Note 
that L is independent of fcmin (though it increases with the size of V{x)). 

The method presented here can be applied to slowly decaying potentials, 
requiring only that V{x) is smooth and [^(a;)! < C/{x)^. This rate of decay is 
the slowest decay rate preserving the shape of I. 

1.1 Heuristics and Intuition 

In this section, we heuristically explain our decay conditions on V{x). In partic- 
ular, we sketch an argument why |V(a;)| < C/{x)'^ allows the interaction region 
to take the form (1.3), but C now varies depending on the potential (rather 
than simply Heisenberg localization). We roughly follow the construction of the 
half- wave parametrix found in [25], and assume V{x) is smooth and decays at 
the rate \V{x)\ < C / {x)'^ . 

To begin, let (/)(a;, k) solve the (classical) Hamilton-Jacobi equation 

\ %4>{x, y, kf + V{x) = lk^ + V{y) (1.4a) 
dxipix, y,k) = k when x = y (l-4b) 
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The solution to (1.4a) is not single valued everywhere, though it is single valued 
on outgoing trajectories. We then let qo{t,x,y,k) solve the classical transport 
equation (valid only where (f){x, y, k) is single valued): 

dtqo = kdxQo + [dxV{x)] dkQo (1.4c) 

qo{x,y,k) = I{x,y,k) (l-4d) 

with I{x, y, k) a smooth function supported in a small neighborhood (a Heisen- 
berg cell) around (.t, y, k) = {xq, xq, ko). Subject to these conditions, it is shown 
in [25, Chapter 4] that: 

e'"'l{x, y,k)- J e4^(^'''.fe)+*('='+^(?'))]go(t, x, y, k)dk 

= O(fc-M^)) = O (1.5) 

The key to making sense of this is the observation that the characteristics 

of the Transport equation (1.4c) are the classical trajectories associated to the 
classical Hamiltonian [1 / 2)p'^ + V {x) . Thus, e*^* moves electrons along classical 
trajectories, ignoring quantum effects. 

It is easily seen that a classical trajectory beginning at (.tq, fco) is outgoing if 
fco/2 > C/{xq)^ > \V{xq)\ and koxo > (the momentum is pointed away from 
the origin). In one dimension, we can simply compute the phase function: 

(l){x,y,k)= I Vfc' + ^y{y) - 2V{x')dx' 

Since k > C/{xq), and \V{y) - V{x')\ < C/{xq)'^ (for different C), we find the 
term under the square root is always positive, and therefore single valued. 

Thus, trajectories originating in the set {(a;o,A;o) : |fco| > C/{xo),xoko > 0} 
are outgoing, and it therefore makes sense to filter waves from this region. 
Provided we consider waves well within this region (i.e. make C sufficiently 
large), quantum effects can not spread the outgoing trajectories too much. 

This argument fails for (a;o,fco) in T. In this region, (l){x,y,k) becomes 
multivalued and non-simple. This is the region we wish to resolve numerically. If 
V{x) decays more quickly than 0{x~^) we make no gain because even though the 
potential is not spread out, low frequency waves are and therefore the interaction 
region still looks like (1.3). 

If V{x) decays more slowly than 0{x^^), the shape of T changes. For in- 
stance, if |V(a;)| < C/{x), then I = {{x,k) : |A:| < C/(a;)i/2}. We believe that 
our method can be extended to such cases, but neglect this to keep things simple 
(see also Section 2.3). 

If k is restricted to the interval [fcmim ^max]? the volume of J behaves like 
0(fcmax log(fc~j?j^)) as fcmin ^ 0. As remarked earlier, covering X with a single 
rectangle (as is done in typical spectral methods), requires solving the problem 
on a region with volume 0(fcmax/fcmin)- This volume is much larger than the 
interaction region, and wastes computation time and space. 
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The algorithm wc construct here involves covering ^fc^^^ femin with phase 
space rectangles of the form [— 2'"L, 2™L] x [— 2~™fcmax, 2~'"/ci„ax]- More pre- 
cisely, on the region [—2~"^L,2"^L], we do calculations on a lattice with lattice 
spacing 27r/(2~'"A'niax) allowing the accurate resolution of frequencies up to 
2~"*fcmax- This allows the computational volume to behave asymptotically like 
the interaction region (up to a constant). See Figure 1 for an illustration. 

The number L is chosen large enough so that waves with frequency greater 
than 2^"'^^A;niax are outgoing by the time they reach x = 2"^L/2. In practice, we 
generally want the kinetic energy to be at least 5—10 times the potential energy 
by the time we filter. For |V(a;)| < Vax~'^, this requires (fcmax/2)^/2 > VqL^^ 
01 L > k~l^^^8Vo/5 (see Figure 2 which illustrates this criteria). This will 
ensure that the phase space filters remove primarily outgoing waves. For the 
interior propagator to work, wc also require that Lfcmax = 0{ln{S^^)) with Si 
a desired error (this ignores logarithmic prefactors; see (2.1c) for the details). 
Last but not least, L must be sufficiently large so that [—L/2, L/2] encompasses 
the region of interest, i.e. the region where we wish to resolve all frequencies. 

Remark 1.1 In spite of many attempts. (1.4) makes a poor basis for a numeri- 
cal method except in the limit k oo (see, e.g. [5]). The reason for this is that 
when (j){x, y, k) becomes multivalued, most techniques break down. As far as the 
authors are aware, no numerical scheme performs better than operator splitting 
spectral methods in regions of phase space with complicated interactions. 



1.2 The Interior Solver 

The interior solver is based on the split step method. In the usual split step 
method, the approximation 



^-i(-{l/2)A+V(x))t _ g«(l/2)AA-t/2 



n 

j=0 



g-*(l/2)AA-t/2 (^]__g^ 



is used. The operator e'ti/^)^"^' is approximated by Fourier transforming the 
wavefunction (using the FFT, or Fast Fourier Transform), multiplying by e^^ /^** 
and inverse Fourier transforming. Of course, the standard Fourier transform is 
based on localizing ip{x,t) in a rectangular region of phase space. 

Our interior propagator is similarly based on (1.6), but the frequency domain 
operators are computed differently. Instead of using a single rectangular region 
of phase space, we cover the interaction region by a union of rectangles (c.f. 
Figure 1). That is, we take a region [— L, L] and use sample spacing 5x to repre- 
sent the region [— L, L] x [— fc^ax, ^max] in phase space. We simultaneously study 
the regions [-2"L, 2"L] x [-2-"fc„iax, ^-"'k^^^] with lattice spacing 2"'6x. To 
deal with the non-uniform sampling in x, we decompose the wavefunction into 
a sum of pieces, each of which is uniformly sampled and localized in a single 
rectangle in phase space. We then use ordinary spectral methods to apply the 
frequency domain operators to these pieces, and reconstruct the wavefunction 
by (spectrally accurate) interpolation afterward. 
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Like the regular FFT, our algorithm achieves spectral convergence. By spec- 
tral convergence, we mean that the all the error comes from parts of the wave- 
function located outside the region of phase space we are considering. We pro- 
vide a rigorous proof of it's accuracy^ in Appendix A. Combining this algorithm 
with standard operator splitting yields a method comparable to normal spec- 
tral methods. The complexity of the algorithm is 0(Mfcinax-^log(fcmax-^)), with 
M = 0(log(fc~j?j^)) the number of dyadic scales needed. 

The dominant part of the error comes from low frequencies, for which is 
too large to resolve even on the coarsest grid (having width 2^2L). However, 
"finite" propagation speed^ combined with virial identities suggest problems 
associated to low frequencies will take an exponentially long time (in the number 
of scales M) to appear; see Section 2.2. 

The approach given in this paper is specialized to the case where the wave 
operator is — (1/2)A, or uj{k) = (1/2)A;^ in 1 space dimension, though general- 
izations are straightforward. Provided operator splitting is valid, the method 
described here can be used. Of course, the rate of downscaling must be sufficient 
to completely capture the interaction region (see Section 2.3). 

1.3 Phase Space Filters for Outgoing Waves 

It should be noted that waves will almost always propagate outside the inter- 
action region; we are just uninterested in them because their behavior is well 
understood. We must therefore come up with a means of removing the outgoing 
waves before they introduce errors into the numerical solution. This is in sharp 
contrast to parabolic equations, where the dynamics of the equation dissipates 
high frequencies (which is the basis for multigrid techniques). 

Since the computational boundary is defined in phase space rather than 
simply position, all the usual methods of open boundaries [2, 3, 11, 12, 16, 28, 6] 
do not apply. Dirichlet-to-Neumann boundary conditions must be imposed at a 
particular curve in x. Artificial dissipation [4, 19] terms like complex potentials 
or the PML will either fail to dissipate high frequency outgoing waves (if placed 
near x = C/fcmin) or will incorrectly dissipate low frequency waves which are 
still interacting with the potential (if placed near x = C/fcmax)- 

The only approach we are aware of which can be generalized to non-rectangular 
regions of phase space is the Time Dependent Phase Space Filter (TDPSF) 
[24]. The TDPSF algorithm (in one dimension, on a rectangular region of phase 
space) is as follows. 

First, we solve (1.1) on the finite region [—L, L]. Inside the regions [— L, —L/2] 
and [L/2,L], we periodically in time apply a phase space filter (with period 
Tstep)- We decompose ip{x,nTstep) = ipout{x) +'iPr{x). The piece ipoutix) is 
strictly outgoing, being localized on the region [L/2,L] x [/cmin, fcmax] in phase 

^Wc neglect floating point errors, and all errors associated to the FFT besides aliasing in 
X and k. These arc small in practice, so wc believe this is reasonable to do. 

^The Schrodinger equation does not have finite propagation speed. However, if frequencies 
are bounded above by fcmax, then velocities are bounded by femax as well. 



5 



Spectral Covering 




Position 



Multiscale Covering 




Position 



Figure 1: A schematic illustration of the phase space regions where different 
algorithms work. The arrows indicate the direction of the flow. 
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Figure 2: The regions of phase space where the filter acts (with 3 scales). The 
arrows indicate the amplitude and direction of motion. The region of phase 
space in which waves interact with the potential is marked. Regions of phase 
space with fc < correspond to incoming waves, and should not be filtered. 



space. The remainder consists of waves located inside the region [—L/2, L/2], 
as well as low frequency waves which may be located in [—L, —L/2] U [L/2, L]. 

Here, "low frequencies" refers to frequencies smaller than A:i„in = 0{{L/2)~^). 
Since our filtering region has width (in position) L/2, we can localize in fre- 
quency no better than 0{{L/2)'^). 

For our problem, with a non-rectangular region of phase space, we wish to 
do the following instead. We will filter off waves in the region [2"^L/2, 2"^L] x 
[2~™fcmax/2, 2~™fcmax] for m = 1 . . . M (for some integer M). Figure 2 diagrams 
the interaction region, filter regions and lattice spacing used. 

1.4 The Main Algorithm 

The general framework of the Multiscale TDPSF propagation algorithm can now 
be described in general terms; the specific details of the internal propagator and 
outgoing wave filter will be treated in Sections 2 and 3 respectively. 
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Algorithm 1 (Multiscale TDPSF Propagation Algorithm) 

Fix an initial condition iJjqIx), and a number of scales M . Set n — 0. 

1. For times t G [nTstcp, {n + l)Tstop]7 solve tl){x,t) on the finite domain 
using operator splitting. The frequency domain operator is approximated 
by Algorithm 5 (see Section 2). 

2. When t = nTgtep, apply phase space filters, replacing il>{x,t) by: 



M 



^OUT 

11=0 



V'(a;,nTstep) 



The projections are calculated using Algorithm 7 (see Section 3). 
3. IfnTstep > Tmax then stop, otherwise goto step 1. 
The computational complexity of this algorithm is 

O(iVMTmaxlogAr) = 0(Tmaxfcr„axMlogfcn,ax). (1-7) 

The parameter N is defined as N = 2L/Sx = 2Lkaia.x/TT. 

Step 2 has this complexity, since we compute M phase space projection op- 
erators, each of which has complexity 0{N/A\ogN) (see Algorithm 7, Section 
3.2). 

It is shown in the Section 2 that Step 1 has complexity 0{NMTstep^ogN). 
Since the loop runs no more than Tmax/^step times, we obtain the correct asymp- 
totic complexity. 

1.5 Connections to other fields 

It is worth mentioning two important areas that our work relates to. First, 
interior solvers similar to ours have been constructed to solve the heat equation 
[15] and the equation V— = / [14]. The major difference between these 
works and ours is that they use the natural decay of the differential operator 
(for large k) to suppress high frequencies, while we use filtering. Those works 
also use the non-uniform FFT, while we make do with the standard one. 

We also remark that similar ideas are used in Coarse Grained Molecular 
Dynamics (CGMD) [21, 9, 8]. In CGMD, one wishes to study the dynamics 
of individual atoms in a small region of space, which are coupled to large scale 
effects (well modclkid by continuum equations) elsewhere. The general ideas are 
similar to ours, although of course our model is continuous at all levels. 

1.6 Disclaimer 

The algorithm presented here is intended to be a proof of principle, and is not 

intended to be directly competitive with other schemes except when fcmax/fcmin 
is large. Additionally, the proof of accuracy of the differential propagator used 
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makes an important assumption that is not strictly true for any real implemen- 
tation. We assume that the only errors made by the FFT are caused by aliasing 
in X and k. We neglect errors caused by floating point arithmetic and errors due 
to discrctizing the integral involved in computing the Fourier transform. Since 
aliasing errors are the dominant part of the error in practice, we believe this 
assumption is reasonable. 

2 Multiscale Calculation of Differential Opera- 
tors 

In this section, we describe a multiscale spectral propagation algorithm (hence- 
forth abbreviated MSP). The MSP uses a non-equispaccd grid to a allow the 
calculation of the Fourier transform of functions which are appropriately local- 
ized in phase space. The crucial fact is that the potential interacts with waves 
only on the phase space regions \k\ < C j \x\. 

The basic idea of the algorithm is the following. First suppose we want 
to apply S'(iV) to a function /(a;), and suppose f{x) has high frequencies lo- 
cated only in the region [— L/2, L/2]. Suppose wc compute the Fourier trans- 
form of x(3^)/(2^)) with x(a^) a smooth partition of unity that is equal to 1 on 
[— L/2, Z//2]. This Fourier Transform is equivalent to Fourier transforming f{x) 
for high frequencies. If we then subtract these high frequencies from /(x), the 
result has only low frequencies. But since the new function has low frequencies 
only, we can Fourier transform it with a larger sample spacing (and reduced 
computational cost) with no loss of accuracy. 

The MSP introduced below uses this by first localizing high frequencies on 
small regions of space, and then by using ordinary spectral methods to propagate 
them. An important technical point is that spectral methods are applied to 
different regions of space with different lattice spacing, and interpolation must 
be used to recover the low frequencies on regions of space with fine sampling. 

2.1 Multiscale Approximation of Differential Operators 

We now describe an algorithm for approximating a differential operator S'(iV) 
applied to smooth functions suitably localized in phase space. In particular, we 
assume that for large x, j{x) has no high frequencies remaining. 

To begin, define the box = i], with L chosen sufficiently large. Let 
^max be a sufficiently large frequency. Both fcmax and L must be large enough 
to satisfy (2.1c). In what follows, 5i is a small parameter, to be adjusted later 
(see Theorem 2.7). 



Xo(a;) = 



-v/ttct 



2 ^_,./,. 




3L/4)]] (2.1a) 
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Po{k) = I [evi{a{k + 3k^^/8)) - erf(a(fc - 3k^^/8))] (2.1b) 

8erfc-^(5i) L ,^ ^ , 

; — < o- < 1 (2.1c) 

Aimax ~ ~ 8erfc"^(5i) 

Here, ★ represents convolution. The standard deviation a is chosen so that 
Xo(a; ^ [-5L/6, 5L/6]) < Si and xo(-4L/6 < x < 4L/6) > 1 - 5i. This means 
that Xo(a^) approximates a partition of unity (supported on [— L/2,L/2]), with 
error 5i. Similarly, Po(fc) approximates a partition of unity, equal to 1 on 

[-fcmax/4, fcmax/4] and smaller than Si on [-fc„iax/2, fc„iax/2]'^. 

Remark 2.1 The function erfc(j;) has the following bounds [1, page 297-298] 
for a; > 0: 

2 e~^^ 2 e~^^ 
—= , < erfc(x) < , (2.2) 

This implies that erfc~^((5) has the bound (for 6 < erfc(l) < 0.158): 



erfc-^(^) < ^/\n{6-^)+\n{^^)/2-ln{2) = 0{^/\n{6-^)) (2.3) 

The functions Xo{x) and Po{k) are the base projections; we now define the 
scaled versions of them: 

Xm{x) = Xo(2-"x) (2.4a) 
P„(fc) = Po{2"'k) (2.4b) 

Before continuing, we formalize our assumption on the phase space localiza- 
tion of f{x). 

Assumption 1 For all m, f{x) satisfies: 

11(1 - Prn{k)).f{x) - Xm{x){l - Pm{k))xm{x)f{x)\\^, < 5^ ||/(x)||^. (2.5) 

One should interpret the operator Xm,(a;)(l — Pm{k))Xm{x) as a "projection" 
onto the phase space region [-2™L,2'"L] x [-2-™fcinax/2, 2-™fci„ax/2]'^. The 
boundaries are of course fuzzy, due to the vagaries of phase space localization. 

Taking the union of these sets for m = O...00, we find that they are a 
covering (using rectangular boxes) of the hyperbola \k\ < C/{x). 

We sample the region Bq with spacing 5x = 27r/3A:niax- This is sufficient to 
resolve frequencies no larger than 3fc,„ax/2 (we explain shortly the reason for the 
extra spacing). We will then sample the regions \ Bm at the rate 2"^~^^Sx. 

This is sufficient to resolve the frequencies 2~'"~^A;max which may exist in the 

region Bm+i, with aliasing errors at most e. 

We add an additional assumption on the function f(x). 

Assumption 2 For all m, S{iW)f{x) satisfies: 

II^(^V)X.„(X)(1 - Pm{k))Xm{x)f{x)\\^.^j,o^^^^ < 62 \\f{x)\\^. (2.6) 

That is to say, S{iV) can not spread waves out too far. 
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This assumption is satisfied by S{iV) = e'^^/^^^"^* for St sufficiently small, pro- 
vided the problem does not have frequencies which are too large. This assump- 
tion is also satisfied even for singular propagators such as ^(iV) = e'l'^l**. 

Remark 2.2 Results very similar to Assumption 2 are proved in detail in 
[23, 26] for S{iW) = e'^^/^'"^^* (in particular, see Section 6). The general 

flavor of the result is as follows. Suppose that the mass of ip{k,t) outside 
the region [— fcmaxj fcmax] is bounded by e. Then for a filter of width w, if 
6t < Cwln((5i)/fcinaxj then S2 is bounded by: 

S2 < C^/Lk^5i + C'e (2.7) 

In fact, the result proved in [23, 26] is more general. That result applies to 
A'' dimensions (replacing VL/cmax by (Lfcmax)''^^), as well as non-differential 
propagators: instead of q^C^/^)^^*^ one can use e~'^^^* with H = — (l/2)A+y(a;) 
for V{x) supported (mostly) inside [—L/2, L/2]. 

We make one further assumption. 

Assumption 3 We assume that \S{k)\ = 1. 

This assumption is not strictly necessary, but it makes the proof simpler. Since 
we are interested in wave propag ators of the form S{k) = e*'^('=)*, this is a 
perfectly reasonable restriction. 

We are now prepared to discuss the algorithm. When discussing the com- 
putational complexity, let N = 2L/5x = 2Lkmax/'^- 

Definition 2.3 The operator Ti,x is the N -point Fast Fourier transform algo- 
rithm operating on the points {-{N/2)5x,-{N/2 + l)5x, {N/2)Sx}. The 
computational cost is 0{N log N). The inverse is denoted by T'^^ ■ 

Algorithm 2 (Phase Space Localization) This algorithm decomposes a 
function f{x) into components which are well localized in phase space. 

1 For m = 0, define 

f+{x) = Xo{x)Ti^\l - Po{k))Ts.Xo{x)f{x) (2.8a) 

f^{x) = f{x)-f+{x) (2.8b) 

The cost of using the FFT region Bo is 0{N log N), while multiplication is 
0{N). No errors are caused by using the FFT with periodic boundaries due 
to the fact that (1 — xo(a;)) vanishes near the boundary of Bq. 

Additionally, define /o^(fc) = for \k\ > fcmax- This is done to control dis- 
cretization errors. 
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2 For m > 1, define: 



(2.9a) 



fmix) = fm-l{x) - .f^{x) 



(2.9b) 



While computing the operator (1 — Pm{k)), we set (1 — Pm{k)) = for \k\ > 
(2/3)2~™fcinax to control discretization errors. 

3 Return the list of functions [foix), fi{x),..., f^ix), fj^ix)]- 

The computational cost of this algorithm is 0{M N log N). This follows 
because for m = . . . M , we compute an FFT of a region with N points at 
a cost 0{N log N). In addition, we must compute an FFT of a region with N/2 
points at a cost 0{N/2logN) ^ 0{N log N). 

RemEirk 2.4 Algorithm 2 will NOT work for an arbitrary function f{x). This 
works for f{x) satisfying Assumption 1 simply because the downsampling used 
on the coarse grids cannot alias high frequencies which are not present. 

Remark 2.5 In steps 1 and 2 of Algorithm 2, we set high frequencies {\k\ > 
2~'"fcmax on each scale) equal to zero. On the level of infinite precision arith- 
metic, this appears unnecessary. In a real computer, floating point errors occur 
which can cause long time instability of the numerical solution. Oversampling 
the grid by a factor of 3/2 and filtering the high frequencies appears to solve this 
problem. This works because floating point errors can be assumed to have arbi- 
trary frequencies, and wc therefore remove (1/3) of them per timestep, whereas 
the effect on the true solution is negligible. 

We have the following result as far as the accuracy of Algorithm 2, which is 
proved in Appendix A.l. 

Theorem 2.6 Suppose that f{x) satisfies Assumption 1. Let f^^'^ he the dis- 
cretization of f^, computed according to Algorithm 2 with M scales. Suppose 
further that a, fcmax and L satisfy (2.1c). Then: 



This result shows the error is linear in the number of scales, and linear 
in fcmax- We conjecture that the dependence on fcmax could be removed; see 
Remark A. 5. 

Before we continue, we provide an interpolation method which provides spec- 
tral accuracy. 




(2.10b) 



(2.10a) 
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Algorithm 3 (Spectral Interpolation) Let f{x) be a funetion on Bm with 
lattice spacing 2™'Sx. f{x) must also be localized in frequency on the region 

\_')—mh, O—mh. 1 
[ ^ "rnax 5 '^maxj • 

Then Tmf{x) is a spectral approximation to f{x) on B^-i with lattice spac- 
ing 2'^~^6x. Calculation ofJmf{x) is done as follows. 

1 Compute the Fast Fourier Transform ofxm-i{x)f{x) on Bm with lattice spac- 
ing 2"^6x. 

2 Compute the inverse Fast Fourier Transform on B„i with lattice spacing 
2™~^Sx. Frequencies with \k\ > 2^'"fcmax are populated with zeros before 
computing this. Discard lattice points outside B^-i- The result isXmf{x). 

Since this algorithm is merely an FFT and an inverse FFT ( on grids with N 
and 2N lattice points, respectively), the computational cost is 0{3N log N) = 
0{N log N). 

We now describe the propagation algorithm. 

Algorithm 4 (Multiscale Spectral Propagator) 

1 Compute [/o^ (x), . . . , f^ix), /M(a;)] by means of Algorithm 2. 

2 For each m, compute S(k)T2^5xfm '^^ '^^ S{k)T2>^6xfM- This has com- 
putational cost 0{M N log N). 

3 Compute the functions: 

9m{x) = ^2'^s^S{k)J^2'^SxfM + •^2M5x'^('^)^2*<f5x/M (2.11a) 
9m-l{x) = J^~J:_i^^S{k)J^2m-iSxfm-l+^m-l9m{x) (2.11b) 

The function gm{x) (defined on Bm) approximates 

M 



S{k) 



S{k)Pmf 



k—m 

on the region Bm with lattice spacing 2™Sx. 

4 Return the function g{x) given by g{x) = gm{x) for x € Bm \ Bm-i- 

We have the following result as to the accuracy of Algorithm 4. 

Theorem 2.7 Suppose S{k) satisfies Assumptions 2 and 3, and f{x) satisfies 
Assum,ption 1. Suppose also that (2.1c) holds. 

Now let g{x) = goix) as calculated by Algorithm 4, but using exact Fourier 
transforms on K instead of the FFT on truncated regions. Let g'^{x) he the 
discrete approximation to g{x) (calculated, with space and frequency truncation, 
but with no integration/ floating point error). 
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Then we have the following bounds: 



Sik)fix) - g\x)L, < \\S{k)f{x) - g{x)\y + \\g{x) - /(a;)|L, (2.12a) 



\\S{k)f{x)-g{x)\\^. 

< Si (M3 + 10M2 + (55/2 + A;„,ax)M + 225A;^a. + 18) ||/||^. 

+ 2||PM/||i2 (2.12b) 



These two formulas should be interpreted as follows. 

Equation (2.12b) is the error caused discarding those parts of the solution 
which Assumption 1 says are small. On the coarsest scale, we discard waves 
which might be significant even if Assumption 1 is true, causing the presence of 
the last term. 

Equation (2.12c) is the error caused by approximating the Fourier transform 
by the discrete FFT algorithm. This is caused both by the fact that our localiza- 
tion procedure introduces tails in x and k (the first term in (2.12c)j and the fact 
that S{k) spreads waves out beyond the limits of truncation (the second term). 
The last term is caused by spatial truncation errors in the lowest frequencies. 

We now present the propagation Algorithm for equation (1.1). 

Algorithm 5 (Multiscale Schrodinger Solver) Take input ipo{x), and fix 
5t > 0. This algorithm approximates il){x,n6t) for n an integer. 

1. Iterate over n = 0..nmax = ^max/^i- 

i Define tpn,i{x) = e^^^*/'^^^^/^^^'t()n{x) . Calculate this using Algorithm 

I 

a Define ^pn,ui^) = e*'**^(^)'(/'„,i(a:;). 

Hi Compute ^n+i{x) = e^^^*/^^^^/^^^tpn,ii{x) , again using Algorithm 4- 

This algorithm has complexity 0{nmaxNM log N). This can be seen since the 
application of Algorithm 4 in steps (i) and (Hi) have complexity 0{NMlogN) 
while step (ii) has complexity 0{NM). 

This algorithm provides 0{5t^) accuracy (0{5t^) for time- dependent poten- 
tials). 

This algorithm is merely the standard 2'nd order operator splitting method. 
The only difference is that wc now use the multiscale propagator to approximate 
the e'(''*/2)(V2)A s^gp instead of the usual FFT. 



with: 



and 



g{x)-g\x)\\^,<5AbM^+ 15 + 




(2.12c) 
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2.2 Interpreting the error bound 

The standard spectral propagator, namely T'^^e^^'^^T^x is said to be spectrally 
convergent for the following reason. Assume f{x) decays exponentially in x and 
/(fc) decays exponentially in k. Assume that e*^*/(x) is localized within [— L, L] 
up to an error of size ^2 ||/|li,2 , where [— i, i] is the computational box. Then the 
spectral propagator can be applied to compute e*'^*/(a;) with error proportional 
to b\ + 82 at computational cost 0{N\ogN), with N = 0{\og{S^^ j). 

Let us discuss why Algorithm 4 provides this level of accuracy. 

First, observe that the computational complexity is 0{MN log N), with N = 
Lk^s,x- Rearranging (2.1c), we obtain the constraint 64(erfc~^((5i))^ < Lkmax- 
Examining (2.3) shows that erfc~^((5i) = 0{^^/log{l/ 61)) for small Si, which 
implies that: 

N = Lk^,^ = 0{\og{6^^)) 

This is spectral convergence m N = Lfcmax; i-e. the error caused by discretizing 
the problem decays exponentially as computational cost increases, provided of 
course that the function is localized sufficiently well in phase space (though not 
exponentially decaying in x and fc, as before). 

There are two additional parameters present: M and \\PMf\\i^2- As M 
increases, we must decrease 5i correspondingly, causing only a logarithmic in- 
crease in N. 82 depends strongly on S{k), and is analogous to errors caused by 
aliasing (waves leaving the right side of the box and entering the left) . 

The last term is a smooth projection onto waves with frequency k < 2~*^A;max- 
This can be bounded if we assume some smoothness of f{k): 



\\PMf\\l2 < Si 



/ 



|fc|>2-'"fe„ 

<sl\ 



J2 



m 



dk- 



1^2 H~ 2 /^max 



-2-Mfc, 

/>) 

-M 



k/2 



max / 2 

2 



dk 



1 1,2 



+ 2-^k^,^C\\{x)f{x)\\l, (2.13) 



The constant C comes from the Sobolev embedding theorem^, since 



m 



Vc 



fik) 



< 



/C|l(x)/(x)||^.. 

For Schrodinger equations, we expect that ||(a;)/(a;)||^2 ~ {p)t (see Sec- 
tion 2.2.1 below). This means that the error due to low velocities behaves 
like 0{2~^/'^kliayi)i which yields exponential decay in M. Therefore, we have 
spectral convergence in N and M; linear increases m N = Lkma,^ decrease 
exponentially, while linear increases in M decrease the error due to low frequen- 
cies exponentially. The remaining error is due to the wave propagator moving 
waves outside the phase space region of interest, which is present in any spectral 

Ul(>tllO(l. 



^Note that more work is required in higher dimensions, since H^(M?'^) L°°(R^'*). 
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2.2.1 Why do we expect ||(a^)/(a;)||^2 to be bounded? 

When used to propagate e*^* as in Algorithm 5, ,f{x) will be an approximation 
to Tp{x,nSt) with ip{x,t) solving the Schrodinger equation. For the Schrodinger 
equation, the virial theorem says that: 

\\{xMx,t)\\l. < \\{xMx,0)\\l^+t'mxMli (2.14) 

In physical terms, this is merely the statement that (x) < (xq) + {p)t, with {p) 
the expected value of momentum. For other types of wave equations the precise 
definition of momentum needs to be changed, but similar identities usually hold. 

Moreover, provided 62 is small, this provides us with a general idea of how 
long a muiicrical simulation (using Algorithm 1) will remain valid. Using (2.14), 
we see that provided y/kmaxUxo) + {p)t) is small, the error will remain 
controlled. This implies that if error e is desired, using 

M = 0(loge-i+log((a;o) + (p)t)) 

scales will control the error. 

In practice, this appears to be an overestimate. Moreover, it suggests that 
as we increase M, wc can use Algorithm 1 for < i < T„iax = 0(2^*^/2). This 
behavior is observed in numerical experiments; see Section 5.2 and Figure 6. 

2.3 Possible Improvement 

Although algorithms 2 and 4 are more efficient than FFT-bascd spectral meth- 
ods for large M, they fail to be competitive for small M (in particular M < 5). 
Careful optimization may reduce this, but made no attempt to do this since our 
main concern is efficiency as M becomes large. 

A more serious problem is that while dyadic downsampling works for V{x) ~ 
Cx~'^ (near x = 00), it will not work for V{x) ~ —Z/ \x\ since the interaction 
region is {{x,k) : \k\ < C/{x)^^^}. For sufficiently large m, Algorithm 1 will 
attempt to filter waves near x — 2™L, k = 2~™A;niax- But the energy of these 
waves is (2~'"A;max)^ — C/(2'"L) < 0, implying that Algorithm 1 will attempt 
to filter waves with negative energy. This is a serious problem, stemming from 
the fact that the interaction region is larger in this case. For these problems, we 
propose using boxes of size Bm = [— 2^'"L, 2'^"^L] each with lattice spacing 2"^5x 
(implying that there are 0(2™) lattice points per scale instead of a constant 
number). Although this is more costly, we believe it will resolve the problem 
for the case of Coulomb potentials (and is still cheaper than using a rectangular 
phase space region). 

3 Phase Space Filtering 

3.1 Introduction to Phase Space Filtering 

In this section we describe the phase space filtering algorithm which will be 
applied. We will apply the phase space filter on a region having size L/A on 
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each side of the grid. This choice is almost certainly not the most efficient, but 
it is the simplest. Consider the following algorithm (see [24] for details): 

Algorithm 6 (Simplified TDPSF Propagation Algorithm) 

Define iul{x) to be zero outside [—L,—L/2], and to he 1 on some region 
inside [— 5Z//6, — 4L/6], andwR{x) — wl{—x). Suppose that ■WL,R{k) is localized 
in the frequency band narrower than [— fcmax/4, fcmax/4] (this is possible provided 
L is large enough). 

Further, define x±(A;) = 1 for k G [±fcmax/4, ±oo), for k < fcmax/8 and 
smooth in between. 

1. For t e [nTstop! (" + l)Tstop]; solve (1.1) by using the Split Step Fourier 
method on [—L,L]. We assume the initial condition is localized in the 
region [-L/2,L/2]. 

2. At times nTgtep, replace ^'(a;, nT^tep) by: 

[1 - WL{x)x-{k)wL{x) - WR{x)x+{k)wR{x)]'^{x,n%tep) 

The operator WL{x)x-{k)wL{x) is a projection onto rightward moving 
waves located to the right of x = L/2, and the wii{x)x+{k)wji{x) is a 
similar operator. 

The time Tsicp is chosen to be (i/4)/fcmax; to guarantee that waves with 
velocity kmax or slower cannot pass through the region [L/2,3L/4] before 
being filtered. 

In practice, we simply take Wl{x) = X[_l+vV7 -L/2-v1^r7] (^) with 
★ denoting convolution. This doesn't quite satisfy the localization condition, 
but comes close enough for practical purposes if e is small enough and L is large 
enough. 

This algorithm is capable of approximating the solution to (1.1) on the 
region [—L/2, L/2], provided ip{x,t) has no outgoing waves with frequency 
k < k^i,^/4:. This can be seen for the following reason. A wave with veloc- 
ity k E [fcniax/4, fcmax] Can travel from [L/2,3L/4] before being removed by the 
application of the operator [1 — WL{x)x-{k)wL{x) — Wii{x)x+{k)]. Since the 
wave never reaches the boundary at a; = L, the boundary conditions are irrele- 
vant, and periodic boundaries can be used. 

This is a simplified version of the Time Dependent Phase Space Filter 
(TDPSF) constructed in [23, 26, 24]. The version in [23, 26, 24] is more ef- 
ficient in many ways. Most notably the buffer region is taken to have width w 
(a parameter independent of L) rather than merely L/2, and can filter waves 
with frequency as low as (Ine/w), with e the desired error. In addition, the 
approach of [23, 26, 24] works in multiple dimensions, while Algorithm 6 does 
not. 

However, the critical problem with both the approach of [24, 23, 26] and 

Algorithm 6 is that they cannot accurately filter low frequency waves. In fact, 
this problem is shared by nearly all methods of open boundaries, including 
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absorbing potentials [19], the PML [4] and Dirichlct-to-Neumann boundaries 
[10, 3, 2, 28, 27, 22, 17, 16]. The reason for this varies depending on the method. 
Dirichlet-to-Neumann boundaries, in all but the simplest (free wave) cases, are 
constructed by means of high frequency approximations, and therefore fail for 
low frequencies. Absorbing potentials and the PML have problems with reflec- 
tion/lack of filtering for small k, analogous to errors made by the TDPSF. 

3.2 The Multiscale TDPSF 

Instead of increasing the width of the filter, thereby lowering /cmin, we simply 
downscale the filter so that it covers the edge of the interaction region. That is, 
we apply the same filtering procedure (with the details clarified), but addition- 
ally using filters with w{2~"^x) and X±(2'"fc) replacing w{x) and X±(fc). With 
M scales this allows us to filter outgoing waves with k as low as 0(2~*^fcmax)- 
We now construct the Multiscale TDPSF. Let w± (x) be a smooth partition 
of [±L/2,±L], and let w±,„{x) = w±(2-"x). Similarly, let X±,n{k) = X±(2"A;). 
We define these operators precisely soon. The operator 

pOUT ^ w+^n{x)x+.n{k)w+^n{x) + ,„ (x)x-,„ (A;)^- ,„ (x) (3.1) 

is then an approximate projection onto outgoing waves. For concreteness, we 
fix a tolerance e, and define: 

w+{x) = (l/2)[erf((a; - (L - 6))/\/2) + erf((a; - {L/2 + b))/V2)] (3.2a) 

b = erf-i(e) = 0{^ln{e-^)) (3.2b) 

b < L/12 (3.2c) 

with w_ (x) defined similarly. The constant b is chosen so that (x) < e to the 
left of L/2 and to the right of L. For this choice to be reasonable, we require 
that L > 2b. The function x+(a;) is given by 

X+{k) = erfc((fc - hn^^/2 - b')/2) (3.3a) 

b' = (3/2) erf-^(e/L) = 0{^/ln{L/e)) (3.3b) 

b' < (3.3c) 

This means that x+{k) > 1 — e when k > 2b'. This cutoff is chosen so that 
even after multiplication by w±{x) (which corresponds to convolution in the k 
domain) PQf{x) will have no more than e frequency content below k = 0. 

The constraints (3.2c) and (3.3c) provide the ultimate limits on the error 
of the phase space filter method. As e decreases, b and b' increase, thereby 
requiring L and &max to increase. The number of lattice points required is 

= 0(Lfc„,ax) > 0{bb') > 0(ln(e))4. 

The operators P,''/'^'"^ can be calculated as follows. 

*Our analysis actually suggests 0{\/ln(e^^)'^ + ln(L) In(e^i)), but the distinction only 
matters when the region of interest is extremely large compared to the accuracy desired, e.g. 
L > . We believe this dependence is an artifact of our analysis: when studying convolutions 
we take absolute values, bounding |2(L/4 - 26) sinc((L/4 - 2&)fe)e-'^^/2| (L/4-26)e-^^/2 
This is suboptimal, but good enough for our purposes. 
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Algorithm 7 (Phase Space Projections) 

Take as input f{x), a function defined on an arbitrarily large domain. This 
algorithm computes P^^'^ f{x). Assume that N lattice points are used to sam- 
ple the region [—L,L\, and that the sampling rate decreases dyadically on later 

regions. 

1. Truncate the domain to [-2"L, -2"L/2] U [2"L/2, 2"L]. 

2. Com,pute w+ji{x)f{x) on [—2"X, —2"L/2] only, and W-^n{x)f{x) on 
[2"L/2,2"i] on/y. 

3. Take the FFT of 

^±,n{x^f{x)j (ind multiply by Xib,n(^)* Note that the 
FFT of each region is done separately. 

4- Invert the FFT and multiply the result by w±^n (x) ■ This yields [Pn^'^ f] {x) 
for X e [-2"i, -2"L/2] U [2"L, 2"i/2]. 

5. For all other x, approximate P^^'^f{x) by 0. 

The computational complexity of this algorithm is 0{N log N), which can be 
seen simply by noting that the FFT has this complexity. All other steps are 
0{N) or 0(1). 

The number of lattice points in Algorithm 7 is only N/A because the region 
[2"L/2,2"L] is sampled with sample spacing 2"(5x rather than simply 5x. 

We arc implicitly assuming that the maximal frequency found in [2"L/2, 2"L] 
is 2~"fcmax rather than fcmax- This assumption actually holds because waves 
with frequency larger than 2~"fcmax are filtered before they can reach the n'th 
layer. 

3.3 Accuracy of the Filtering Algorithm 

In [23, 26], it is proven rigorously that if V{x) « in the filter region. Al- 
gorithm 6 is accurate^. We show that if the solution 'ip{x,t) has frequencies 
boimdcd between K and fcmax and if V{x) is well localized, a phase space filter 
of width 0{\n{e)/K) will provide accuracy of order O(Tniaxe) (ignoring loga- 
rithmic prefactors). In [24, 23, 26], K is taken to be fcmin, while in this work 
we take K = fcmax/2 (leaving low frequencies to be filtered by the downscaled 
filters). These error bounds are proven only for rapidly decaying V{x). Proving 
accuracy for slowly decaying V{x) is more difficult, and involves showing that 
Pn^'^ projects accurately onto outgoing waves of —(1/2) A -|- V{x) (a nontrivial 
problem of scattering theory). 

To demonstrate the accuracy of a single phase space filter, we solved (1.1) 
using Algorithm 6. The initial condition was compu- 
tational region was x G [—51.2,51.2] and k G [—31.4,31.4] (1024 lattice points), 
and Tjnax was taken to he AL/k (enough time for waves to wrap around the box 

^We actually prove the result for a variant of Algorithm 6 which extends to multiple 
dimensions. 
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Phase Space Filter Errors 




Frequency ($k$) 

Figure 3: A graph of error (||V'(a;, t) - V-dCa;, i)llL2([-LX]) / llV'o(a;)|ji2) vs fre- 
quency of an outgoing pulse for the phase space filtering algorithm. 

twice). A phase space filter with varying e, but fixed L and fcmax, was used to 
remove outgoing waves. Waves with frequency fc > 13 are accurately filtered, 
while low frequency waves are not. The results are plotted in Figure 3. 

4 Justification of the Method 

We discuss briefly why Algorithm 1 works. 

Begin by assuming that V'(a;,0) = iphix) + ^pB,+ {x) + ^pB,~{x), where tphix) 
is localized on [— L/2, L/2], ijjB^+ix) is localized on [L/2, L] and has frequencies 
primarily k > fcmax/2, and ipB.-ix) is localized on [L/2, L] and has frequencies 
primarily on fc < A;niax/2. (For simplicity, we consider waves on the right side 
only, the left being treated identically.) 

The component tphix) can be propagated safely for at least a time Tstcp = 
L/4fcinax- This follows since the maximal velocity is fc,„ax, and it would take a 
time 2Tstep for to cross the buffer region. 

The function ipB,-{x) can be propagated safely as well, at least for time 
L/(fcmax/2) — STstop- This is the time it takes for waves with velocity fcniax/2 
to cross a region of width L. The region in question is the low frequency buffer 
that is added to the edge of the grid (the second scale) . 

Finally, ipB,+ {x) can not be propagated safely, since it consists of high fre- 
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qucncy waves which arc about to enter the second scale; however, the first filter 
removes tpB,+{x) before it enters the second scale and corrupts the solution. 

The same argument can be repeated for arbitrarily many scales, until 2~™fcinax < 
/Cmin- Thus we see that if fcmin is known, we can solve the problem using 
m = 0(log(A;max/fcniin)) scales. 

4.1 What if A;inin is unknown? 

In the event that fcmin is unknown, we can still solve the problem if Tmax is known 
and finite. Suppose we have M scales. Then outgoing waves with frequencies in 
the range [2~^k^ax, ^max] can be successfully filtered. The only waves which 
can cause a problem are waves with k G [—2~'^k^a^^,2~^kjaa.^]- Supposing 
the initial condition to be localized on [—L/2, L/2], these waves can travel no 
farther in space than L/2 + 2~ ^ kynaxTmax- Due to their low frequency, these 
waves are located on the coarsest scale, which extends as far as x = 2^"^ L. 

This implies that the slow waves hit the boundary when 2~*^/CmaxImax ~ 
2"^L. We therefore choose m so that 

2^^L 

T^nax < -T (4.1) 

SO that the low frequency waves do not have sufficient time to reach the bound- 
ary. Thus, with an unknown fcmin but known (large) Tmax, the choice M ^ 
log(TniaxfcmaxT~^) guarantees the accuracy of the simulation. The complexity 
of the algorithm is therefore 

0[{T^ax/St) l0g(T,„ax)fcmaxi log(fcmaxi)] (4.2) 

It should be noted that some Dirichk;t-to-Neumann implementations for 
the Schrodinger equation achieve time-complexity 0(Ti„ax log Tmax) [16, 17, 18] 
as well (for compactly supported potentials). Very roughly, where Dirichlet- 
to-Neumann boundaries require storing the time-history on the boundary, we 
require storing extra low frequency outgoing waves. 

5 Numerical Results 

In this section we present the results of numerically implementing Algorithm 1. 
The algorithm was implemented in the Pjrthon programming language, using 
the libraries Numarray, Matplotlib and FFTW [13]. The code is available at 
http://cims.nyu.edu/~stucchio under the Gnu Public License. 

5.1 The Free Schrodinger Equation 

We begin by running some comparisons to solutions of the free Schrodinger 
equation (namely (1.1) with V{x,t) = 0). We solve this equation taking the 
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Error vs Frequency, 3 Scales 




5 10 15 20 



Frequency 



Figure 4: A graph of error vs the velocity of an outgoing pulse. ^()e{x,t) is the 
exact solution, ip{x, t) the numerical one. 

initial condition tp{x, 0) = (47r(7)^/2gifexg-xV2<T^^ ^^^^ fc = 1, 2, . . . 21. With this 
initial condition, the exact solution is 

exp (ifc(x - fct/2)) / -\x-kt\'^ \ 

M^^t) - ^1/42^1/2(1 + ,i/^2)i/2 1,2^2(1 + zt/a2) ) ■ ^^-^^ 

We solve (1.1a) by means of Algorithm 1 using 3 scales. Each scale has 1024 
lattice points, taking Sx = 0.1 on the finest scale. The parameter e was taken to 
be 10~®, leading to b' — 7.88 (as per (3.3)). The maximal resolvable frequency 
on the finest scale is k = 41.8. Outgoing waves are filtered from the n'th scale 
at fc = 2-"13.06. The timestep is taken to be 5t = 2"^. 

In each simulation, Tmax is taken to bo Tmax = 204. 8/A;, which is more than 
enough time for the outgoing wave to reach x = 51.2 and return to the origin. 
The quantity 

E{k)= sup ||V'(a;,t) - V'e(a;,t)||i2r_25.6,25.6l (5-2) 

te[0,T„ax] 

was computed, and the result is plotted in Figure 4. 

Figure 4 shows that the error remains uniformly below 10~^, for a simulation 
with fj = 4. The errors will get progressively worse with high frequencies, 
however, but this can be resolved by increased the rate of sampling. 

It should be noted that for /c « 1, if we ran the simulation out to time 
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Figure 5: A graph of error vs the width of an stationary pulse. ipeiXjt) is the 
exact solution, tjj{x,t) the numerical one. 



t ~ 400+, we would see errors due to these waves wrapping around the compu- 
tational domain. This problem is resolved by the use of additional scales. 

Another set of simulations was run, this time varying a with (7 = 1,..., 128. 
The maximal time was Tmax = 50.0 in all cases. Note that this initial condition, 
when (T = 2^ = 64 is not even localized inside the region [—51.2,51.2], and is 
therefore completely inaccessible by standard absorbing boundary techniques on 
a box of size [—25.6, 25.6] since the typical wavelength of this solution is longer 
than the computational domain. The results are plotted in Figure 5. 

For this simulation, wc could have solved the problem with cr = 128 by using 
a grid with 1024 lattice points and dx = 0.4. The reason for using the multiscale 
TDPSF algorithm is that low and high frequencies can be solved simultaneously. 
Indeed, tests involving initial conditions containing both wide gaussians (with 
large a) and fast ones (with k S [1, 15]) achieved accuracy of 10~^ as well. This 
is important for problems where low and high frequencies mix. 



5.2 The Long Range Schrodinger Equation 

In atomic physics one often considers potentials V{x,t) may be decaying very 
slowly in x. The prototypical case is V{x) = —Z/ \x\, though we restrict our- 
selves to the simpler case V{x) x~'^ near ±oo: 

V{x) = — . + 20e--'/^ 

1 + 25.6-2 Ixj^ 
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The initial data is taken to be (47r)~^/^e"^ 1'^''^ . The fine grid is taken to occupy 
the region [-102.4, 102.4] with the first filter occupying the region [-102.4, -51.2]U 
[51.2, 102.4]. The simulation was run with 2, 3 and 4 scales up to a time t = 1500. 
To generate an "exact" solution to compare to, the same simulation was run us- 
ing the standard FFT on the large region [-52428.8, 52428.8] (requiring 524288 
lattice points). The results are plotted in Figure 6. 

To give a picture of the dynamics of the solution, iIj{x, t) (for various values 
of t) and V{x) are plotted in Figure 7. The parameters are 5x = 0.1, dt = 
and 6 = 10~^. On the finest scale, gaussians with |6fco| > 9.1 are filtered, 
while others are left on the coarser scales. As can be seen from Figure 7, the 
solution is rough on the interval [—50,50], and high frequencies are present. 
If fine sampling were not used, aliasing errors would occur. The waves which 
escape the potential have very low velocity {v w 2) when they reach x = 50, too 
slow to filter by normal methods. 

As can be seen from Figure 6, more scales allow the simulation to run for a 
longer time. The times at which the simulations begin to fail arc slightly before 
t = 200 (with 2 scales), t = 400 (with 3 scales) and t = 800 (with 4 scales). This 
is in rough agreement with (2.13) and the discussion in Section 2.2, although 
the agreement is not exact due to the presence of the potential ((2.14) is valid 
as written only when V{x) = 0). The two-scale error occurs more quickly than 
would be expected since waves move rapidly at the bottom of the potential (by 
direct analogy to a classical particle). 

6 Conclusion 

In this work wc have presented an algorithm for approximating the solution of 
time-dependent dispersive wave equations on by domain truncation. All 
other methods we are aware of attempt to solve the problem on a rectangular 
region of phase space, which is inaccurate for waves with large wavelengths. 
We argue that a hyperbolic region {{x,k) : \k\ < C/{x)} in phase space is 
the correct region to consider, and provide a spectrally accurate algorithm for 
approximating differential operators on this region. 

Additionally, we extend the Time Dependent Phase Space Filter algorithm 
[24] to accurately filter outgoing waves regardless of frequency. The computa- 
tional complexity is proportional to log(fc~j'jj) rather than k^l^ unlike the PML 
or absorbing potentials. 

It should be noted for compactly supported potentials in up to two space 
dimensions, Dirichlet-to-Ncumann boimdaries do obtain [16, 17, 18] similar 
0(Ti„ax logTmax) Complexity (recall Eq. (4.2)). We are aware of no compa- 
rable results for long range potentials, however. 
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Errors vs Time, multiple scales 
lo-if ' ' ' ' ' 




Time 

Figure 6: A graph of error (relative error, measured in L^) vs time for a test 
involving a long range potential. As predicted in section 4.1, the time that the 
simulation remains accurate increases with the number of scales. 
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Figure 7: Plots of |V'(a;,t)| for various times, as well as V{x). For t > 40, the 
motion remains confined to the bottom of the potential well. 
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6.1 Nonlocal Wave Operators and Other Spectral Prob- 
lems 

The method described here can be apphed to other equations, inclu ding tho se for 
which the wave operator is nonlocal. The wave operator uj{k) = ^/k^ + |J? — ij is 
an example take from relativistic quantum mechanics (with // > 0) and quantum 
field theory (with = 0). Since \/—A -\- jj? — ijl is nonlocal, neither finite differ- 
ences nor finite elements can be used to solve such equations. There is nothing 
precluding the use of the MTDPSF for these purposes, and prelimi nary num eri- 
cal tests suggest Algorithm 1 works for propagation when u!{k) = \/k'^ + jj? — jjL. 

Additionally, the techniques of Section 2 can be used for other spectral prob- 
lems. For instance, we have solved the Poisson equation on R using Algorithm 
5 with S{k) = — = k"^ (using the quadrature rule from [7] to deal with the 
singularity at fc = 0) as well as the equation y/—Au = /. While a new Poisson 
solver is unremarkable, it is notable that extension to non-local equations is so 
simple; usually clever tricks are required [14]. 

6.2 The Missing High Frequencies 

The Multiscale TDPSF algorithm approximates il^{x,t) only on the interior of 
Bq. On the interior of B„i, the munerical solution approximates Prnijj{x, t). This 
means that on the vast majority of the computational domain, the information 
that the algorithm yields is incomplete. 

In the broader sense, however, we believe that the TDPSF provides all the 
relevant information for the problem. The philosophical assumption that under- 
lies open boundaries is the assumption that outside a fixed region of space, the 
behavior of the solution is free and uninteresting. However, as discussed earlier, 
this assumption is conceptually flawed with regard to low frequencies. Low fre- 
quencies interact over a proportional to their wavelength, which is potentially 
much greater than the computational box (see Figure 2). 

For this reason, we argue that the MTDPSF algorithm is providing the 
right information, while all other algorithms are conceptually flawed when low 
frequencies are present. 

6.3 Future Directions 

When the MTDPSF removes waves from the computational grid, they are re- 
moved because they are outgoing and moving nearly freely. In some applications, 
knowledge of their motion after they leave the domain may still be desired. 

This can be rectified at low cost, however, and we intend to investigate this 
in future works. In the high frequency limit, solutions to (1.1) have simple 
behavior. A WKB expansion in the time- variable yields a sequence of transport 
equations which can be safely truncated for frequencies which are high relative 
to the size of the potential (this corresponds to the absence of turning points). 

By the time the TDPSF has removed high frequency waves, they have moved 
into the WKB regime. This means that if further information concerning their 
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propagation is desired, it can be obtained cheaply by using WKB and then 
solving a transport equation. Thus, we believe the MTDPSF algorithm can 
form a building block for solving (1.1) on extremely large regions of space. 



A Proof of Correctness 

In this section our ultimate goal is to prove the correctness of Algorithm 5. 

We first prove the accuracy of Algorithm 2 in A.l. We then use the results 
there to prove the accuracy of Algorithm 4. 

Our main tool is a result describing the accuracy of computing a differential 

operator by means of a Fourier transform on a finite box. Essentially, the result 
states that if we restrict a function to a finite box —[L,L]^, then the error in 
approximating S{iV)(p{x) is the mass of S{iV)(p{x) located outside [—L,L]^. 
In addition, if we apply a frequency cutoff (such as the one caused by sampling), 
then there is an additional error corresponding to the mass of S{iV)(j){x) which 
is cut off. 



Theorem A.l Let S{iV)ip{x) satisfy the hypothesis of the Poisson summation 
formula, that is \S{iS/)ip{x)\ < C{x}^+^ and S{ik)(f{k) < C(fc)'^+^ LetS{k), 

Sh{k) be continuous bounded Fourier multiplication operators which are equal 
for k £ B (where B is some closed set). 
Then: 



S{iV)ip{x) - e'''''-^/^Si,{nk/L)<f>{TTk/L) 



< ||5(iV)<p(f+2Ln)||jy,((j_^_^]«)C)+ Hk) 



sup 



S{k) - Sb{k) 

(A.1) 



Remark A.2 A result which is similar to this one appeared in [20], where it is 
used to show that systems on a large enough box can exhibit transient radiative 
behavior (over short times) in the same way that systems on can. 



Proof. The Poisson summation formula states that: 

^ f{x + n2L)= e'''^-^/^f{-Kk/L) 



(A.2) 



We let /(A?) = S{k)(p{k). Then, by rearranging (A.2), we find: 
S{iV)<fi{x)- Y e'''''-^/^S{nk/L)<fi{nk/L) = - ^ S{i^)(p{x + 2Ln) (A.3) 
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Since S{k) and Sb{k) arc equal on B, wc can add and subtract 
E,rfc/i,ez« e^''^-^/^{Sb{-Kk/L) - S{TTk/L))ip{Trk/L) to both sides of (A.3), to ob- 
tain: 



(A.3) 

= - S{iV)ip{x + 2Ln) + e'^'^-^/^I^Sbi-Kk/L) - S{Trk/L))(p{wk/L) 

= ^ (S'6(iV)-S'(iV))(/?(f +2Ln) (A.4) 

where the last line follows by applying (A. 2). We now take norms and apply 
the triangle inequality. We find that: 

Y \\S{iVMx + 2Ln)\\H. = \\S{iW)^{x + 2Ln)\\Hs((^_L,L]^^a^ 
and that: 



Y ll(^6(«V) - S{'i^))^{x + 2Ln)\\jjs = \\{Sb{iV) - S{iV))ip{x + 2Ln)\\jjs 



< 



m 



sup 



S{k) - Sb{k) 



Wc put everything together to obtain the result we seek. 



□ 



Using this result, we will show that appropriate differential operators can 
be reasonably calculated. The main assumption which must be satisfied is that 
S'(iV) can not move Pm{k)xrn{x)f{x) outside Bm- If S{iV) is a local operator, 
this property is immediately satisfied. 

As a warm-up, we prove the accuracy of Algorithm 3. 



Theorem A.3 Let f{x) be a function such that: 



III 



\k\>2- 



Xm-lf{k) 



dk < 



The 



\\1mf{x)-f{x)\\^.^^^_^^<{e + 5,)\\f\\^. 
Recall 5\ is a desired error hound, and Xk depends on it (c.f. (2.1)). 



(A.5) 
(A.6) 



Proof. Apply Theorem A.l with S{iV) = 1 to the box Bm-i- In this case, 
note that Sb{k) = 1 for |fc| < 2~™fcinax, and otherwise. This shows that: 



\\S{k)Xm-lf - Sb{k)Xm~lf\\ L^{B^_,) 



<e 
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The result of this is a function supported on [— 2~'"~'^fcniax, 2~'"~'^fcmax]- Af- 
ter padding 2~"'~^fcmax < \k\ < 2~'"A:„iax with zeros, the result is a function 
supported on [— 2~™fcmax, 2~™fcmax]- Inverse Fourier transforming yields a func- 
tion g for which \\g — Xrn-i.f\\]^2(^g < £• By the definition of Xm-i (provided 
(2.1c) is satisfied), we find, that |^y7i_i - 1| < (5i for X G Bm- Thus: 

\\f - St,{k)Xm-lfh^B^) 

< 11/ - S{k)Xm-lf\\mB^^ + \\S{k)Xm-lf - St{k)Xm-lf\\mB^) 

< 11(1 - Xm-i)fh^B^_,) + e ll/IL. < {6, + e) 11/11^. 
This is what we wanted to show. □ 



A.l Correctness of Algorithm 2 



In this section it is our goal to show that Algorithm 2 correctly approximates 
the list of functions [foik), /j^(fc), . . . , /;j^(^)) /m(^)]- Recall the definition of 
Xk and Pfc and a from (2.1). 



Lemma A. 4 Suppose that 



m 



, < '5i||/(x)||^.. Then the 

i ( [ — fcmax , fcmax] ^ ) 

function (1 — Po)xof = fo can be accurately approximated on the box [—L,L] 
with maximal frequency Amax- 



\{l-Po)xf-{^-Po)xf\ 



LHl-L,L]) 



< ^1 2 + 



2lT(T 



11,2 



(A.7) 



Here, Pq is the operator Pq restricted to the box [—L,L], with a frequency cutoff 

at /Cmax • 

Proof. By the aliasing theorem, for k < kj^ax, the discrete version of Pq 
agrees with the continuous version of Pq, whereas for k > k^ax, they differ by 
2. Theorem A.l implies therefore that: 

||(1-P0)X/-(1-P0)'^X/|L.([-L.L]) 

< Wxf - x/IIl2([_l,l]) + \\Poxf - Poxf\\m[-L,L]) 

< + \\Poxf\\m[-L,L]0) + 2 \\xf\\m[-k^,^,k^,^]o,dk) 

< \\PoXf\\m[-LM^) + 2^1 Wmh^ (A.8) 

The last inequality follows by assumption. Treating the first term requires 
somewhat more work. 

The operator Pq can be written (in the x-domain) as a convolution, 



Pn = 



an 



1/2 



^ 2A;iiiax sinc(fcn;iax^) 
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where a satisfies (2.1c). Thus, we need to bound: 



/ dx 


[ dx' 


l[-L,L]C 


Jr 



^ ^ ^ 2fcniax Sinc(fcinax(2; X )) 



Xo{x')f{x') 



/dx I dx' 



-{x—x')'^ /a 



xo{x')f{x': 



< 



4\-'a-'ki,fdx'xo{x'f\f{x')f I dx 

Jr J[-l,l]c 



<<a. 11/(^)112= 



Xo{x') 



dx- 



1/2 



-{x-x'f/a^ 



(A.9) 



The second Une follows from the first by noting that |sinc(2;)| < 1 for z real. 
The function inside the L°° norm can be explicitly calculated: 

^1/2^1/2^ [erf(a-i(3L/4 + x')) + evf{a-'{x' - 3L/4))] 

(erfc(2i/2^-i(L - a;')) + erfc(2i/2^-i(a;' + L))) 

The quantity inside the square brackets is 1xq{x), and is therefore bounded by 
2^1 for x' ^ [— 5iy/6, 5L/6]. For x' E [— 5L/6, 5L/6], the quantity inside the 
erfc function is at least as large as y'2L/6a < L/3a; this implies (if cr, L, fcmax 
are chosen according to (2.1c)) that the quantity inside the round brackets is 
bounded by 2Si. Thus: 



oi,2 

iA.9)<d,^\\f{x)\\l. 



Substituting this bound on (A.9) into (A. 8) yields the result we seek. 



(A.IO) 



□ 



Remark A. 5 We believe the factor of fcmax present in (A. 7) could probably be 
removed by a more careful analysis. However, the cost of our laziness will only 
be a logarithmic factor at the end of the day. 



Proposition A. 6 We have the hound: 



16fcn 



2-K(J 



11/(^)11 



1,2 



(A.11) 



Proof. Suppose that ^/^(a;) — fm''^{x)\ij2 < -Em- Then simply note that: 



m+l 



r-,d 
J m+l 



|[1 - Xra{x){l - P™)Xm(a;)]/™(x) - [1 - Xm{x){l - P?k)Xm{x)]f+^''{x)\\^, 
< II ([1 - Xm(x)(l - Pm)Xm{x)] " [1 " Xm{x){l - Pi)xM]) /-(x)||^, 

+ ||[1 - xm{x){l - P:^)xm{xm-''{x) - f-{x))L^ (A.12) 
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Since Xm{x){l — P:^)xm{x) is a self-adjoint, positive operator bounded by 1, we 
find that: 

II 1 - Xm{x){l - -P™)Xm(a^) 11^(^2,^2) < 1 

Thus: 



(A.12) <6,(2+ 1 ||/(a;)||^2 + ||/-'''(x) - f-{x)\\^. 



By summing over the Em, (and summing the geometric series 2 from m = 
. . . oo to avoid messy dependence on m), we obtain the result we seek. □ 

Proposition A. 7 We have the bound: 

\K{x)-fy{x)\\^.<S, (^2m+i^) ||/(x)||^2 (A.13) 

Proof. A simple calculation, entirely similar to the previous 



Jm+l Jm+1 



L2 



— ||Xm(^)(l Pm)Xm{,X)fm Xm(2')(l P^Xm{x)f^' ||j^2 
< II [Xm(a;)(l - Pm)Xm{x) - Xm{x){l - P^)Xm{x)] /-||^, 



+ ||Xm(a;)(l - Pi)Xm{x) {f-^\x) - f-{x)) I 



i2 



<5, (2m + 2+i^) ||/(a=)||,2 



with Em given as in the proof of Proposition A. 6. □ 

We have now shown that the effect of discretization on Algorithm 2 is min- 
imal, and Theorem 2.6 is proven. 

A. 2 Correctness of Algorithm 4: Proof of (2.12c) 

We prove here the accuracy of Algorithm 4. 

Proposition A.8 (Algorithm 4, Step 2) Let /i+ = S{k)f+ and let h+;'^ 

be the discrete approximation to Define hj^ similarly. Assume also that 
\S{k)\ < 1. Then the following error bound holds: 



hi - h^W^, < 5i ( 2m + ) ||/(x)||^2 + 62 ||/(x)||^2 (A.14) 



16k 
/27rcr 

Note that this result computes the error associated to Algorithm 4 up step 2. 
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Proof. Observe that (with S'^{k) the discrete approximation to S{k)): 

W^m ~ ||i2 = \\^ik)fm~^ {^)fm ||l2 

< \\S{k){f+- fy)\\^, + \\[S{k) - S''{k)]f+\\^, (A.15) 

Proposition A. 7 allows us to bound the first term in the second line of (A.15). 
The second term is bounded by applying Theorem A.l, and noting that As- 
sumption 2 bounds the mass localized outside B^- Thus, we have: 

(A.15) < 5, {2m + ||/(x)|L. + 5, \\f{x)\\^. 



This is what we wanted to show. □ 

Step 3 of Algorithm 4 is exact. Since gfy-i{x) is already band-limited by 
by discretization, spectral interpolation is exact. This implies immediately the 
following result: 

M 

/ I \ 



\gt{x) - gML^^s^) < E ^1 ('^+ ^) 



+ 52 \\f{x)\\j^. + ^'(fc)/M'' - S{k)f]^{x) 



< 5i {2M{M - m) + (M - m)i^) ||/(a;)||^. 

+ 52{M - m) \\f{x)\\^. + 2 ||/m(x)||^. (A.16) 



1,2 



To simplify, we assume that m = (maximizing the right of (A.16)), ob- 
taining: 

(A.16)<5i (^2M2 + M^^) ||/(:r)||^. 

+ ^2M||/(a;)||^.+2||/^(a;)||^, (A.17) 

The term here which is not controlled by a factor oi 61^2 is ||./m(-^)||^2- This 
term corresponds to waves on the coarsest scale, and can be thought of roughly 
&sPM.f{x). The reason is as follows. By Assumption 1, XTO(a;)J-*m(^)Xrn(a;)/(x) ft 
Pmf{x)- This means that each time we subtract Xm(a^)(l — Pm{x))xm{x), we 
are approximating Pm{x); the remainder after this is done M times consists 
solely of low frequencies. 

Proposition A. 9 We have the following hound: 

II/- - P„/„_i||^, <^m + 2)U+ 5, ll/IL. (A.18) 
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Wc will prove this in Appendix B.l (the proof is straightforward but long). 
The basic idea is that high frequencies of f{x) are localized where Xm{x) ~ 1, 
so the effect of Xm{x) is negligible. 

This immediately leads to the following result: 



Proposition A. 10 We have the bounds: 

m 

f;n-I[Pkf 



k=0 



- , 3m^ + 15m 
< Si ( + 6 + 75k^^ 



L2 



1L2 (A.19a) 



m— 1 



k=0 



< 61 (Sm^ + 15m + 150kn,^ + 12) \\f\\^, (A.19b) 



Proof. First, note that: 



Pm ■ ■ - PlPof = Pm ■ ■ - Pi {Pof — /o ) + Pm ■ ■ ■ Plfo 

= Pi(Po/ - /o") + Pm... P2{Plfo - fl) +Pm... ^2/1" 



E 

j=0 



Bringing to the left, taking norms and using Proposition A. 9 (as well as the 
fact that ||-Pm||£(i:,2,L2) ^ 1) ^^'^ simple geometric series bounds^) shows: 



\Pm ■ ■ ■ PlPof - /m||j:,2 < E 



3=0 



L2 



<E0' + 2) 3 



4-2^ 



^1 



, 3m(m + l) + 12(m+l) , 15A;„,ax(i + 2) 

<'^i I o +E 



3=0 



A -23 



im? + 15m + 12 



<5l 



3m^ + 15m 



+ 6 + 75fc„ 



max \\J 11^2 



(A.20) 



^Namely the fact that XlTLo 2"^ < 1/(1 - 1/2) and also ETLo J^"^ < 1/(1 - 1/2)^ 
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This is (A.19a). Finally, note that: 



k=0 



\\fm-l ~ fm~ [Pm-l ■ ■ ■ PlPof - Pm ■ ■ • -Pl-Po/) || 

< II/+ - Pm-l . . . PlPof\\i^2 + \\-f-+Pm... PlPofW^. (A.21) 

Using (A.20) to bound (A.21) yields (A.19b). □ 
We now add and subtract 111=0 P^f inside the last norm of (A. 17), yielding: 



(A.17)<5i f2M2 + M^5^ 
+ <52M||/(a;)||^.+2 



m 



k=0 



k=0 



(A.22) 



We apply Proposition A. 10, in particular (A. 19a) with m = M to (A.22). We 
also observe that HfeLo Pk is a product of commuting, positive operators, and 
therefore WUZoPkfW < II^m/||. This yields: 



(A.22) < di 2M2 + M 



^ , \\f{x)h2+62M\\fix)h.+ 

+ ^1 (3M2 + 15M + 12 + 150fc„,ax) ||/||i2 + 2 \\Pm!\\l^ (A.23) 

We now simplify (A.23) and formalize it as Proposition A.ll. 

Proposition A.ll (Algorithm 4, Step 3) We have the following error 
bound: 



\9m{x) -9m{x) 



< Si IbM^ + 



15 + 



16fcn 



2iTa 



M + 12 + 150fc„ 



+ <52M||/(a;)||^.+2||PM/IL2 (A.24) 
In particular, setting m = in (A.24) recovers (2.12c). 

A. 3 Correctness of Algorithm 4: Proof of (2.12b) 

We have now controlled the discretization errors associated to the calculation of 
g'^{x), and Theorem 2.7 is half proved. The only thing remaining is to estimate 
the difference between g^ix) and S{idx)f{x) on B^, when computed using the 
exact operators. 
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First, observe that we can write: 



fix) = (1 - Po)/ + Pof = (1 - Po)f + (1 - Pi)Po/ + PiPof 

= (1 - Po)/ + (1 - Pl)Pof + (1 - P2)PlPo/ + P2P1P0/ 



= (1 - Po)f + 
Similarly, we can write: 

9o{x) = 



M /m-1 \ 

J2(^-Pm)il[Pkf] 
m=l \/s=0 / 



M 



[]Pfe/ (A.25) 



/s=0 



M 



_m— 



M 



(A.26) 



Wc wish to compute a bound on S{k)f{x) — gaix). We will first estimate how 
close E^^, /+ is to (1 - Po)/+ [Ei=i(l - Pm) (nr=o' ^fc/)] • Using (A.19b) 
term by term, we find: 



M 



m=0 



M 



/' m— 1 



.m=l 



^fe=0 



1,2 



M 



< E '^1 (3™^ + 15to + 150fcmax + 12) 



It2 



< Si{M^ + 9M^ + (20 + 150fc,„ax)M 

+ 150fc„ax + 12fc^ax + 12(M + 1)) 11/11^2 (A.27) 

Combining this with (A.25) and (A.26), and using the fact that S{k) is 
unitary shows that: 



\\S{k)f{x)-go{x)\\^, 



< 



M 



/ m—1 



M 



(i-Po)/+E(i-^-) n^'^/ -E/- 



m— 1 
M 



^fe=0 



m=0 



1,2 



+ 



n^'fe/-5(fc)/M 



fe=0 



L2 



< 5i(M3 + + (20 + 150fcn,ax)M + 150fc„,ax + 12fc„,ax + 12(M + 1)) 11/11^2 

+ \\PMfh2 + \\S{k)f]^\\^, (A.28) 

Using (A. 19a) again, adding and subtracting nm=o inside the last norm 
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of (A. 28), we obtain the bound: 

(A.28) < 5i(M3+9M2+(20+150/e„,ax)M+150fcmax+12/e„,ax+12(M+l)) ||/||^. 



+ \\PMfh2 + <5i + + 6 + 75fc^a.) II/IL. 



M 



k=0 



< 5i{M^ + 9M^ + (20 + 150fc„,ax)M + 150fc„,ax + 12fc„,ax + 12(M + 1)) ||/||^. 

, /3M2 + 15M \ „ „ 

+ 2 ||Pm/IL2 + 5i { + 6 + 75A:^axj 11/11^2 

< 5i (M3 + lOM^ + (55/2 + k^^)M + 22bk^^^ + 18) + 2 WPmIWl^ 

(A.29) 

Thus, we have proved (2.12b) as well. 

B Technical Points 

B.l Proof of Proposition A. 9 

We are now almost ready to prove Proposition A. 9. We first state a Lemma 
used in the proof. 

Lemma B.l Provided (2.1c) holds, and m < j are integers, we have the fol- 
lowing bound: 

||[(1 - Pj{k)) - Xj{x){l - Pj{k))Xj{x)]Xm{x)Pm{k)Xm{x)\\^(^.^^.^ 

The idea behind Lemma B.l (proved shortly) is that [(1 — Pj{k)) — Xj(x)(l — Pj{k))xj{x) 
is "almost" supported on S^, while Xm{x)Pm{k)xm{x) is supported on Bm- 
Since m < j, Bm C Bj and so Bm n Bj^ is empty. Of course, this is only 
approximate., so it must be quantified. 

Proof of Proposition A.9. Note that/" = [l-Xm{x){l-Pm{k))Xm{x)]f~+i. 
We wish to show that (A. 19a) holds. Begin by writing: 

- P^f^-l = [1 - Xm{x){l - Pm{k))Xmix)]f~_i - Pmfm-1 
= [1 - P„(fc) - Xm{x){l - Pm{k))Xm{x)]f-_, 

= [l-P™(fc)-Xm(a;)(l-P„(fc))xm(a;)](l-Xm-i(a;)(l-P™-i(fc))xm-i(a;))/-_2 

= [1 - P„(fc) - Xm{x)il - Pm{k))Xm{x)]f-^2 

- [1 - Pmik) - Xm{x){l - Pmik))Xmix)]Xrn-l{x)^ f~_2 

+ [1 - Pm{k) - Xm(a;)(l - Pm{k))Xm{x)]Xm-l{x)Pm-l{k)Xm-l{x)fm-2 

(B.2) 
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Define rm-j{x) by: 

Tm-j = -[1 - Pm{k) - Xm{x){l - Pm{k))Xmix)]Xm-j{xf f~-j-i 

+ [1 - Pm{k) - Xm(a;)(l - Pm{k))Xmix)]Xm-j{x)Pm-j{k)Xm-j{x)f~_j_i 

If we then expand /^_2) and repeat ad-nauseum, (B.2) becomes: 

(B.2) = [1 - P„(fc) - Xm{x)il - Pmik))Xm{x)]f-_2 + Tm-l 

= [l-P„(fc)-Xm(a;)(l-P™(fc))xm(a;)](l-Xm-2(a;)(l-P„-2(fc))xm-2(a;))/-_3 

+ rm-1 + rm-2 

m 

= [1 - Pm{k) - Xm{x){l - Pm{k))Xm{x)]f{x) + ^ Tm-j (B.3) 

3=0 

The term f{x) surfaced since f^ix) = (1 — Xo{x)Poik)xoix))f{x). Assume for 
the moment (it will be proved shortly) that: 

||r„-i(x)||^. < (3 + (15/c^ax/4)2— ) 6, ||/||^. . (B.4) 

Then: 

m 

||(B.3)|L. < ||[1 - P„(fc) - xm{x){l - P^{k))xUx)]f{x)\\^.+J2\\rm-jh2 

3=0 

<5,\\f\\,. + {m+l) (3 + ^^) 6,\\f\\^. 

<(m + 2) (3+^)^1 11/11^. (B.5) 

The term - Pmik) - Xm{x){l - ^^m(fc))Xm(a;)]/(a;)||^2 was bounded merely 
by applying Assumption 1. Thus (A. 18) is proved, pending a proof of (B.4). 

Boundedness of Vm-j- Proof of (B.4) 

We wish to compute a bound on: 

rm-3 = -[1 - Pni{k) - Xm{x)(l - Prn{k))Xrn{x)]Xm-3 {xf .f~_j_i 

+ [1 - Pm{k) - Xm{x){l - Prn{k))Xm{x)]Xm-jix)Pra-j{k)Xm-j{x)f~_j_.^ 

By Lemma B.l, we can bound the last term on the right of (B.2) by 

(1 + fl^) ^1 \\f^-3-A^ < (1 + ^) ^1 Wfh^ ■ 

The second to last term is bounded as follows. The function Xm-jix)"^ f^_j_i 
has L^.norm smaller than ||/||^2 outside D = [-2™-J5L/6, 2'"-^'5i/6], since 
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Xm{xY < there. Inside the region [— 2™4L/6, 2'"4L/6], we can write Xm{x) = 
1 + h{x), where \h{x)\ < 6l Since D C [-2™4L/6, 2^4^/6], we need only con- 
sider this region. Here, we find: 

||[1 - Pm{k) - Xm{x){l - Pm{k))Xm{x)]Xrn-jixy.f:rn-j-l\\L2^D) 

= ||[1 - Pm{k) - (1 + h{x))il - P™(fc))(l + /l(x))]x™-,(x)V~-,-l||^.(^) 



< [hix){l - Pm{k))il + h{x))]Xm-j{xrf, 



-i-l|lL2(D) 



+ \\[il + Hx))il-PUk))h{x)]Xrn-j{x)'f-_^_,\\^,^^^ 

< 2 \\h{x)\\^^ \\Xm{x){l - Pm{kmc^r.^L^) ||Xm-i(x)V;;-,-i|L.(^) 

< 2<5i ll/IL. 

This completes the proof. □ 

Before proving Lemma B.l, we require one additional result to be stated. 
Lemma B.2 Let g{x) = (21 a)er'-^' /(^'o^) ^ let m < j be integers. Then: 

Xm {x 

Proof of Lemma B.2. Note that: 

i'^-Xjix)) < ^-2HL/6-2HL/6]c{x) + Sili_2j4L/6-2UL/6](x) (B.7a) 
Xmix) < l[-2"5L/6,-2'"5L/6] + '5ll[-2'"5L/6,-2™5L/6]^' (B.7b) 

Note that the distance between [-2HL/6, -2HL/6f and [-2"5i/6, -2'"5L/6] 
is at least 2J(4L/6- 2™--'5i:/6) > 2^(4^/6 -5L/12) = 2^1/4 (since m-j < 1). 

We now break g{x) = gi{x) + g2{x), with .91(2;) = g{x) for \x\ < 2^ L/8 and 
zero otherwise, and g2{x) = g(x) — gi{x). Simple integration shows that: 

\\92{x)\\^, < 2erfc < 2erfc(L/8a) < 25^ (B.8) 

This follows since a > L/8erfc~^(Ji) and erfc(a;) is monotonically decreasing in 
X. By Young's inequality, this implies that ||s'2(a;)i»r||£(£,2 ^,2) < 25i. 

Now, letting f(x) = (1 — Xj{x))g{x) * Xm{x)h{x) (for some h{x) G L^), we 
observe that: 

\f{x)\ = Ji^- Xj{x))g{x - x')^Xm{x')h{x')dx' 

(1 - Xj{x)) [go{x - x') + giix - ,/)] * Xm{x')h{x')dx' (B.9) 
Substituting the bounds (B.7) into (B.9), we find: 

(B.9) < [l[_2J4L/6,-234L/6]c(a;) + (5il[_2J4L/6,-2J4L/6](a;)] 

X (ffo(a:)*+3i(a;)*) 

X [l[-2'"5L/6,-2'"5L/6](a;) + <5ll[-2'"5L/6,-2'"5L/6]c(a;)] /l(a;) (B.IO) 
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Note that (B.9) is a produet of the form {Ai + SiBi){A2 + (5ii32)(^3 + (Si-Bs). 
When the sum is expanded, there will be 8 terms, only one of which does not 
have a factor of 6i in it. This term is 

l[-2J4L/6 -2ML/6]c(a;) [gaix) * 1[_2™5L/6,-2™5L/6] (a^)'i(a;)] 

But since the support of goix) has width 2^L/A, which is the distance between 
the support of the first characteristic function and the second, this term must 
be zero. 

There are another three terms of order 5i (one of which is actually 2Si , the 
term coming from gi{x)-k), and 4 more terms of order Sf. Provided < 1/4, 
we find that the sum of these terms is 5^i. This is what we wanted to show. □ 

Proof of Lemma B.l. A calculation. First: 



[(1 - Pj{k)) - Xj{x){l - Pj{k))Xj{x)] Xm{x)Pm{k)Xm{x) 
= (1 - X3{x))Pj{k)Xm{x)Pm{k)Xm{x) 

+ Xj{x){l - Pj{k)){l - Xj{x))Xm{x)Pm{k)Xm{x) 

Recalling (2.1a), we observe that |(1 — X3{x))Xm{x)\ < Si. This implies that: 

||Xj(a;)(l - Pj(/e))(l - Xjix))Xrnix)Pmik)Xmix)\\^f^^2^L^j < Si 

Thus we only need to bound (1 — Xj{x))Pj{k)xm{x)Prn{k)Xrn{x). This can be 
done by examining the integral form of (1 — Xj{x))Pj{k)xmix) (the operator 
Pm{}^)Xm{x) has norm bounded by 1): 



\il~xdx))PAk)Xmix)i-)\ 

(1 - x,(x))e-(--')V(2^'^)^^ sine (^^(^ _ 

J{1- Xj{x))e-^''-'''^'/^'''^^\mix')i-)ix')dx' (B.ll) 



/ 



4 •2-' 



The last line follows since |sinc(^r)| < 1. By Lemma B.2, we find that the norm 
of the integral operator is bounded by 5Si, thus: 

'ik 

(BAl)<^5Si\\{.)\\^. (B.12) 
This is what we wanted to show. □ 
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